//////////////
// 1. setup //
//////////////

// Aggregate
var A nu RealR px Pi u R
// Low type
cn_L W_L cu_L delta_L lambda_u_L u_lag_L ut_L u_L theta_L v_L Vj_L W_L mu_L Vv_L Vn_L Vu_H WNB_L
// High type
cn_H W_H cu_H delta_H lambda_u_H u_lag_H ut_H u_H theta_H v_H Vj_H W_H mu_H Vv_H Vn_H Vu_H WNB_H

varexo epsilon_A epsilon_nu;

///////////////////
// 2. parameters //
///////////////////

parameters rho_A rho_nu std_A std_nu vartheta sigma beta M_L M_H alpha epsilon_p phi Pi_ss A_L A_H  p_L p_H Upsilon psi eta mass_L R_ss 
rho_R delta_pi rho_nu std_nu  rho_A std_A; 

// household    
        mass_L=0.5; // mass of L agents
        beta = 0.96^(1/12); // discount factor
        sigma = 2.0 # CRRA coefficient
        vartheta = 0.90; // home production / replacement rate
        alpha = 0.6; // elasticity of MF
        gamma = np.nan; // worker bargaining weight, determined endogenously
        
        //intermediary goods firms
        kappa = 0.58; //vacancy flow cost
        psi = 40; //separation elasticity
        p = np.nan; // maximum separation probability, determined endogenously
        Upsilon = np.nan; // Vj at maximum separations, determined endogenously
        eta = 0.5; // exponential weight on nash wage in wage rule
       //final goods firms
        epsilon_p = 6.0; //price elasticity      
        phi = 1000.0; //Rotemberg cost

       //monetary policy
        rho_nu = 0.0; //persistence 
        std_nu = 0.0001; //standard deviation
        R_ss = 1/beta;
        rho_R = 0.96; //inertia
        delta_pi = 1.5; //inflation aggressive
        Pi_ss = 1;
         rho_A = 0.0; // persistence 
         std_nu = 0.0001; //standard deviation

   // to be determined from data


    alpha 
///////////////////////////////////////
// 3. initial guess for steady state //
///////////////////////////////////////



//////////////
// 4. model //
//////////////

model;

  //////////////////
  // a. exogenous //
  //////////////////

// put shock process here

exp(A)= exp(A(-1))^rho_A*exp(epsilon_A) ;

exp(nu)=exp(nu(-1))^rho_nu*exp(epsilon_nu);

  ////////////////
  // b. workers //
  ////////////////

    [name = 'worker consumption L'] exp(cn_L) = exp(W_L);
    [name = 'worker consumption H'] exp(cn_H) = exp(W_H);

    [name = 'unemp consumption L'] exp(cu_L) = vartheta*W_L;
    [name = 'unemp consumption H'] exp(cu_H) = vartheta*W_H;

    [name='euler L'] (exp(cn_L))^(-sigma) = 
        beta*exp(RealR)*(
          exp(delta_L(+1))*(1-exp(lambda_u_L(+1)))*(exp(cu_L(+1)))^(-sigma) + 
          (1-exp(delta_L(+1))*(1-exp(lambda_u_L(+1))))*(exp(cn_L(+1)))^(-sigma)
        );

  ///////////////////////////////////////////
  // c. labor market dynamics and matching //
  ///////////////////////////////////////////
    
    [name = 'u_lag_L'] exp(u_lag_L) = exp(u_L(-1));
    [name = 'ut_L'] exp(ut_L) = exp(u_lag_L) + exp(delta_L)*(1-exp(u_lag_L));       
    [name = 'tightness_L'] exp(theta_L) = exp(v_L)/exp(ut_L); 
    [name = 'job_finding_L'] exp(lambda_u_L) = M_L*exp(theta_L)^(1-alpha);  
    [name = 'job_filling_L'] exp(lambda_v_L) = M_H*exp(theta_L)^(-alpha);     
    [name = 'u_L'] exp(u_L) = (1-exp(lambda_u_L))*exp(ut_L);      

    [name = 'u_lag_H'] exp(u_lag_H) = exp(u_H(-1));
    [name = 'ut_H'] exp(ut_H) = exp(u_lag_H) + exp(delta_H)*(1-exp(u_lag_H));       
    [name = 'tightness_H'] exp(theta_H) = exp(v_H)/exp(ut_H); 
    [name = 'job_finding_H'] exp(lambda_u_H) = M*exp(theta_H)^(1-alpha);  
    [name = 'job_filling_H'] exp(lambda_v_H) = M*exp(theta_H)^(-alpha);     
    [name = 'u_H'] exp(u_H) = (1-exp(lambda_u_H))*exp(ut_H);


  ////////////////////////////
  // d. final and wholesale //
  ////////////////////////////
  
  [name='Philips'] 1 - epsilon_p + epsilon_p*exp(px) = 
      phi*(exp(Pi)-Pi_ss)*exp(Pi) - 
      beta*phi*(exp(Pi(+1))-Pi_ss)*exp(Pi(+1))*(1-exp(u(+1)))/(1-exp(u));

  /////////////////////
  // e. intermediate //
  /////////////////////

    [name='Bellman for Vj L'] exp(Vj_L) = exp(px)*exp(A)*A_L-exp(WNB_L) + (1-exp(delta_L(+1)))*beta*exp(Vj_L(+1)) - beta*exp(mu_L(+1));
    [name='Bellman for Vv L'] exp(Vv_L) = -kappa + exp(lambda_v_L)*exp(Vj_L);
    [name = 'delta L'] exp(delta_L) = p_L*(exp(Vj_L)/Upsilon)^(-psi);
    [name = 'mu L'] exp(mu_L) = p_L*psi/(psi-1)*Upsilon * (1-(exp(Vj_L)/Upsilon)^(1-psi)) / (1-p_L*(exp(Vj_L)/Upsilon)^(-psi));

    [name='Bellman for Vj H'] exp(Vj_H) = exp(px)*exp(A)*A_H-exp(WNB_H) + (1-exp(delta_H(+1)))*beta*exp(Vj_H(+1)) - beta*exp(mu_H(+1));
    [name='Bellman for Vv H'] exp(Vv) = -kappa + exp(lambda_v_H)*exp(Vj_H) ;
    [name = 'delta H'] exp(delta_H) = p_H*(exp(Vj_H)/Upsilon)^(-psi);
    [name = 'mu H'] exp(mu_H) = p_H*psi/(psi-1)*Upsilon * (1-(exp(Vj_H)/Upsilon)^(1-psi)) / (1-p_H*(exp(Vj_H)/Upsilon)^(-psi));

  ///////////////////
  // f. Worker values //
  ///////////////////
    
   [name = 'Value of employment L'] exp(Vn_L)= (exp(WNB_L)^(1-sigma)-1)/(1-sigma)+beta*exp(delta_L(+1))*(1-exp(lambda_u_L(+1)))*exp(Vu_L(+1))+beta*(1-exp(delta_L(+1))*(1-exp(lambda_u_L(+1))))*exp(Vn_L(+1));
   [name = 'Value of employment L'] exp(Vu_L)= ((vartheta*exp(WNB_L))^(1-sigma)-1)/(1-sigma) +beta*exp(lambda_u_L(+1))*exp(Vn_L(+1))+beta*(1-exp(lambda_u_L(+1)))*exp(Vu_L(+1));
   [name = 'Value of employment H'] exp(Vn_H)= (exp(WNB_H)^(1-sigma)-1)/(1-sigma)+beta*exp(delta_H(+1))*(1-exp(lambda_u_H(+1)))*exp(Vu_H(+1))+beta*(1-exp(delta_H(+1))*(1-exp(lambda_u_H(+1))))*exp(Vn_H(+1));
   [name = 'Value of employment H'] exp(Vu_H)= ((vartheta*exp(WNB_H))^(1-sigma)-1)/(1-sigma) +beta*exp(lambda_u_H(+1))*exp(Vn_H(+1))+beta*(1-exp(lambda_u_H(+1)))*exp(Vu_H(+1));

  ///////////////////
  // g. Wages //
  ///////////////////   

    [name = 'Nash Wage L'] gamma*(exp(Vn_L)-exp(Vu_L))=(1-gamma)*(exp(Vj_L)-exp(Vv_L));
    [name = 'Real Wage L'] exp(W_L) = W_ss_L * (exp(WNB_L)/W_ss_L)^(eta);

    [name = 'Nash Wage H'] gamma*(exp(Vn_H)-exp(Vu_H))=(1-gamma)*(exp(Vj_H)-exp(Vv_H));
    [name = 'Real Wage H'] exp(W_H) = W_ss_H * (exp(WNB_H)/W_ss_H)^(eta);

  ///////////////////
  // g. Output and unemployment //
  ///////////////////

    exp(X)=exp(A)*(mass_L*a_L*(1-u_L)+(1-mass_L)*A_H*(1-u_H));
    
    exp(u)=1-(mass_L*(1-u_L)-(1-mass_L)*(1-u_H));
  ////////////////////
  // h. government //
  ////////////////////

    [name = 'Taylor'] exp(R) = R_ss*(exp(R(-1))/R_ss)^(rho_R)*(exp(Pi)/Pi_ss    )^(delta_pi*(1-rho_R))*exp(nu);
    [name = 'RealR'] exp(RealR) = exp(R)/exp(Pi(+1));

end;

/////////////////////
// 5. steady state //
/////////////////////

initval;
end;

resid;
steady;
check;

///////////////
// 6. shocks //
///////////////

@#include shocks_path